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REPORT  DOCUMENTATION  PAGE  (SF298) 
(Continuation  Sheet) 

Continuation  for  Block  13 


ARO  Report  Number  56225.5-EG 

Pressure  effects  on  the  relaxation  of  an  excited  ... 


Block  13:  Supplementary  Note 

©2015  .  Published  in  The  Journal  of  Chemical  Physics,  Vol.  Ed.  0  142,  (1)  (2015),  (,  (1).  DoD  Components  reserve  a  royalty- 
free,  nonexclusive  and  irrevocable  right  to  reproduce,  publish,  or  otherwise  use  the  work  for  Federal  purposes,  and  to  authroize 
others  to  do  so  (DODGARS  §32.36).  The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and 
should  not  be  construed  as  an  official  Department  of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other 
documentation. 


Approved  for  public  release;  distribution  is  unlimited. 


THE  JOURNAL  OF  CHEMICAL  PHYSICS  142,  014303  (2015) 


* 


Pressure  effects  on  the  relaxation  of  an  excited  nitromethane  molecule 
in  an  argon  bath 

Luis  A.  Rivera-Rivera,1  Albert  F.  Wagner,2  Thomas  D.  Sewell,1  and  Donald  L.  Thompson1 

1 Department  of  Chemistry,  University  of  Missouri-Columbia,  Columbia,  Missouri  65211-7600,  USA 
2Argonne  National  Laboratory,  Chemical  Sciences  and  Engineering  Division,  Argonne,  Illinois  60439,  USA 

(Received  30  September  2014;  accepted  4  December  2014;  published  online  5  January  2015) 


Classical  molecular  dynamics  simulations  were  performed  to  study  the  relaxation  of  nitromethane 
in  an  Ar  bath  (of  1000  atoms)  at  300  K  and  pressures  10,  50,  75,  100,  125,  150,  300,  and  400  atm. 
The  molecule  was  instantaneously  excited  by  statistically  distributing  50  kcal/mol  among  the  internal 
degrees  of  freedom.  At  each  pressure,  1000  trajectories  were  integrated  for  1000  ps,  except  for  10 
atm,  for  which  the  integration  time  was  5000  ps.  The  computed  ensemble-averaged  rotational  energy 
decay  is  ~  1 00  times  faster  than  the  vibrational  energy  decay.  Both  rotational  and  vibrational  decay 
curves  can  be  satisfactorily  fit  with  the  Lendvay-Schatz  function,  which  involves  two  parameters:  one 
for  the  initial  rate  and  one  for  the  curvature  of  the  decay  curve.  The  decay  curves  for  all  pressures 
exhibit  positive  curvature  implying  the  rate  slows  as  the  molecule  loses  energy.  The  initial  rotational 
relaxation  rate  is  directly  proportional  to  density  over  the  interval  of  simulated  densities,  but  the 
initial  vibrational  relaxation  rate  decreases  with  increasing  density  relative  to  the  extrapolation  of  the 
limiting  low-pressure  proportionality  to  density.  The  initial  vibrational  relaxation  rate  and  curvature 
are  fit  as  functions  of  density.  For  the  initial  vibrational  relaxation  rate,  the  functional  form  of  the 
fit  arises  from  a  combinatorial  model  for  the  frequency  of  nitromethane  “simultaneously”  colliding 
with  multiple  Ar  atoms.  Roll-off  of  the  initial  rate  from  its  low-density  extrapolation  occurs  because 
the  cross  section  for  collision  events  with  L  Ar  atoms  increases  with  L  more  slowly  than  L  times 
the  cross  section  for  collision  events  with  one  Ar  atom.  The  resulting  density-dependent  functions  of 
the  initial  rate  and  curvature  represent,  reasonably  well,  all  the  vibrational  decay  curves  except  at  the 
lowest  density  for  which  the  functions  overestimate  the  rate  of  decay.  The  decay  over  all  gas  phase 
densities  is  predicted  by  extrapolating  the  fits  to  condensed-phase  densities.  ©  2015  AIP  Publishing 
LLC.  [http://dx.doi.org/10. 1063/1 .49043 14] 


I.  INTRODUCTION 

Advanced  internal  combustion  engine  designs1  for  efficient 
and  clean  combustion  use  dilute  fuel  mixtures  at  pressures 
considerably  higher  than  currently  employed.  This  strategy 
reduces  peak  combustion  temperatures,  leading  to  lower 
pollutant  formation  while  retaining  power  density.  Pressures 
as  high  as  500  bars  are  within  the  conceptual  range  of  practical 
devices,  far  higher  than  typical  peak  operating  conditions 
of  100  bars  in  current-generation  engines.  A  rational  design 
process  requires  knowledge  of  relaxation  and  reaction  rates  at 
higher  pressures  than  are  typically  studied  experimentally  or 
theoretically.  In  particular,  high  pressure  affects  the  unimolec- 
ular  and  recombination  rates  and  branching  ratios  because 
relaxation  and  reaction  processes  compete  to  control  the  fate 
of  highly  excited  molecule  or  radical. 

Competition  between  relaxation  and  reaction  occurs  in 
many  complex  systems  that  do  not  involve  combustion,  such 
as  the  atmosphere.  That  competition  is  usually  represented  in 
models  by  the  master  equation  (ME).  There  are  a  number  of 
reviews2^4  of  ME  applications  to  modeling  combustion  and 
atmospheric  processes,  and  they  generally  emphasize  that  ME 
modeling  is  not  reliably  predictive  stemming  mainly  from 
the  lack  of  knowledge  of  energy  transfer,  angular  momentum 
effects,  and  potential  energy  surfaces  (PESs).  This  lack  of 


knowledge  concerns  not  only  the  absence  of  specific  relaxation 
data  on  the  many  different,  and  often  transient,  species  in  com¬ 
plex  systems  but  also  the  fact  that  the  implementation  of  ME 
formulations  in  current  models  involves  additional  approxima¬ 
tions.5  More  rigorous  implementations  of  the  ME  approach  are 
an  active  area  of  research6  that  will  grow  more  challenging 
for  higher  pressures  for  which  gas-molecule  collisions  are  no 
longer  isolated  from  one  another.  In  order  to  provide  species- 
specific  relaxation  data  and  to  contribute  to  improved  ME 
formulations,  studies  of  collisional  energy  transfer  over  a  wide 
interval  of  pressures  are  needed. 

A  number  of  techniques  have  been  used  to  measure 
energy  relaxation  in  highly  excited  molecules.  These  include 
infrared  (IR)  fluorescence  detection,7  high-resolution  IR  laser 
transient  absorption  spectroscopy,8  time-resolved  diode  laser 
spectroscopy,9  time-resolved  ultraviolet  absorption  spectros¬ 
copy,10  time-sliced  velocity  map  ion  imaging  in  molecular 
beams,11  and  laser-Schlieren12  measurements.  Most  of  the 
experimental  measurements  using  these  techniques  have  been 
carried  out  at  low  pressures  for  which  a  collision  of  the 
excited  molecule  with  the  buffer  gas  is  rare  and  well  iso¬ 
lated  from  any  subsequent  collisions  of  the  excited  molecule. 
However,  several  studies  have  been  performed  for  higher 
pressures  with  the  objective  of  relating  the  results  to  relaxation 
studies  in  liquids. 13  Of  particular  interest  are  the  high-pressure 
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experimental  studies8161,8^’14,15  that  have  used  supercritical 
fluids  for  which  the  density  can  be  conveniently  varied  from 
small  values  typical  of  lower-pressure  environments  to  large 
values  that  are  only  achieved  at  very  high  pressures. 

Theoretical  studies  of  collisional  energy  transfer  involv¬ 
ing  highly  excited  molecules  have  been  dominated  by  classical 
trajectory  simulations.  These  are  an  alternative  to  rigorous 
quantum  dynamics  calculations,  which  are  not  feasible  for 
systems  of  this  size.  Although  the  neglect  of  interference  ef¬ 
fects  in  classical  mechanics  affects  its  validity  at  low  ener¬ 
gies,  quantum  effects  should  be  of  minor  importance  for  en¬ 
ergy  transfer  at  the  high  energies  of  interest  here  where  the 
quantum-classical  correspondence  principle  should  be  valid. 
The  trajectory  simulations  can  be  logically  divided  into  three 
groups.  In  the  first  group  are  simulations  of  ensembles  of 
isolated  collisions  of  bath  gas  species  with  a  molecule  or 
radical  excited  to  a  specific  internal  state.  Simulation  for  a 
variety  of  initial  states  of  the  molecule  produces  results  that 
can  be  fit  to  the  probability  transition  matrix  required  by  the 
ME  formalism;  examples  can  be  found  in  Refs.  5  and  16-18. 
In  the  second  group,  the  successive-collisions  approximation 
of  Lendvay  and  Schatz19  is  used  in  which  IVR  is  assumed  to 
randomize  all  but  the  conserved  properties  (total  energy  and 
angular  momentum)  in  the  post-collision  state  of  a  molecule 
prior  to  the  next  collision.  No  ME  is  involved,  rather  ensemble 
averages  over  independently  generated  sequences  of  succes¬ 
sive  collisions  are  used.  At  high  enough  pressures,  the  collision 
frequency  will  exceed  the  IVR  rate  and  this  approach  will 
break  down.  In  the  third  group,  molecular  dynamics  (MD) 
simulations  explicitly  follow  the  motion  of  the  excited  mole¬ 
cule  and  all  the  bath  gas  species  as  they  collide  with  each  other 
and  ultimately  relax  to  a  thermal  distribution.  This  is  by  far 
the  most  expensive  of  the  three  groups  of  simulations  because 
the  vast  majority  of  collisions  are  between  bath  gas  species 
that  serve  only  to  sustain  the  bath  gas  thermal  population. 
However,  this  approach  is  the  only  one  that  is  equally  valid 
at  low  and  high  pressures  (up  to  and  including  the  condensed 
phases  of  the  bath).  Because  of  the  expense,  such  simulations 
are  unusual,  especially  for  highly  excited  molecules  of  greatest 
interest  in  combustion  for  which  the  time  required  for  thermal- 
ization  can  be  quite  long.  Several  examples  can  be  found  in  the 
work  of  Heidelbach  et  al. 211  for  azulene/CCL  and  of  Paul  et  al. 23 
for  the  CgFe/Ni  system. 

While  it  is  impossible  to  briefly  summarize  the  vast  body 
of  experimental  and  theoretical  results  on  collisional  energy 
transfer,  there  are  two  general  observations  of  particular  rele¬ 
vance  to  this  paper.  First,  to  the  limits  of  experimental  and 
theoretical  resolutions,  relaxation  of  initially  highly  excited 
molecules  to  thermal  equilibrium  is  typically  characterized  by 
single-exponential  decay.  However,  there  are  experimental10,12 
and  theoretical19,22,23  examples  in  the  literature  where  multiple 
relaxation  rates  have  been  observed.  Dove  et  al.10  measured 
the  collisional  energy  transfer  of  highly  vibrationally  excited 
CS2  in  27  bath  gases  at  300  K.  The  energy-loss  profiles  are  not 
exponential,  including  those  for  all  the  noble  gases.  Lendvay 
and  Schatz19  calculated  the  collisional  energy  relaxation  of 
highly  excited  CS2  in  He  and  Xe  at  300,  1000,  and  2000 
K.  The  relaxation  of  the  molecule  was  found  to  be  exponen¬ 
tial  in  He  but  non-exponential  in  Xe.  Kiefer  et  al.12  reported 


shock-tube  laser-Schlieren  measurements  of  the  high-tempera- 
ture  pyrolysis  of  ethane  and  CF3CH3  in  which  they  observed 
two  vibrational  relaxation  rates  at  temperatures  above  and 
below  the  onset  of  dissociation.  Li  and  Thompson22  calculated 
the  vibrational  relaxation  of  L  in  Xe  at  ~300  K.  At  high  initial 
energy  and  for  an  anharmonic  I2,  the  vibrational  energy  relax¬ 
ation  was  found  to  be  non-exponential.  Paul  et  al.23  calculated 
the  relaxation  of  vibrationally  excited  CgFg  in  N2  bath  at  298  K. 
For  all  pressures  considered  (35  atm-707  atm),  the  relaxation 
of  the  molecule  was  found  to  be  non-exponential. 

The  second  general  observation  is  that  while  measure¬ 
ments  and  calculations  of  relaxation  rates  are  directly  propor¬ 
tional  to  pressure  at  pressures  low  enough  for  isolated  gas- 
molecule  collisions,  the  measurements  of  Schwarzer  et  al. 14 
and  the  calculations  of  Heidelbach  et  al.1'1-1  and  Paul  et  al.21' 
show  that  at  high  enough  pressures,  the  relaxation  rate  grows 
more  slowly  than  in  direct  proportion  to  the  pressure.  The 
Schwarzer  et  al.  measurements  of  the  relaxation  of  azulene 
in  a  variety  of  supercritical  fluids  including  CO2  show  that 
a  change  in  relaxation  rate  growth  from  direct  proportional 
dependence  to  slower  growth  occurs  at  a  density  of  ~  1  mol/L 
or  ~25  atm  (estimated  using  the  ideal  gas  law).  A  companion 
MD  calculation  by  Heidelbach  et  al.2{>  for  azulene/CCL  qual¬ 
itatively  reproduced  both  the  measured  low-  and  high-density 
relaxation  rates,  thus  reproducing  the  departure  of  the  rate  from 
proportionality  to  pressure  at  high  density.  Finally,  the  MD 
simulations  of  Paul  et  al.23  were  carried  out  for  a  variation  in 
pressure  of  over  a  factor  of  20  and  the  computed  relaxation 
rates  show  a  less-than-proportional  growth  of  the  relaxation 
rate  with  pressure  at  high  pressures.24 

We  report  here  the  results  of  MD  simulations  of  the 
relaxation  of  50  kcal/mol  of  initial  excitation  in  nitromethane 
(CH3NO2)  in  a  300  K  Ar  bath  at  eight  different  pressures  from 
10  atm  to  400  atm.  In  these  simulations,  the  molecule  and 
all  the  Ar  atoms  are  explicitly  followed  for  at  least  1000  ps, 
by  which  time  the  nitromethane  has  nearly  reached  thermal 
equilibrium.  The  50  kcal/mol  excitation  is  ~10  kcal/mol  less 
than  the  weakest  bond  of  the  molecule,  namely,  the  C-N 
linkage  between  the  CH3  and  NO2  groups.25  Consequently,  no 
dissociation  occurs  in  the  simulations.  The  pressure  range  of 
10  atm^MX)  atm  straddles  the  gas/supercritical-fluid  transition. 
In  the  phase  diagram  of  Ar,  the  critical  point  is  at  150.68  K  and 
48.0  atm.26  Consequently  at  300  K,  only  pressures  below  48 
atm  are  in  the  gas  phase.  From  48  atm  to  400  atm,  Ar  is  a 
supercritical  fluid.  At  more  than  an  order  of  magnitude  higher 
pressure,  Ar  exits  the  supercritical  phase  and  condenses  into  a 
solid. 

This  is  the  first  report  among  a  planned  series  of  reports 
on  simulations  designed  to  examine  the  relaxation  of  diatomic 
molecules,  polyatomic  molecules,  and  radical  species  in  an  Ar 
bath  up  to  quite  high  pressures.  The  rest  of  this  paper  is  orga¬ 
nized  as  follows:  Sec.  II  contains  a  discussion  of  computational 
methods  including  the  PES  used  in  the  simulation.  Section  III 
contains  a  description  of  the  results  for  both  vibrational  and 
rotational  relaxations.  The  focus  of  Sec.  IV  is  the  discussion  of 
a  model  that  rationalizes  the  observed  pressure  dependence  of 
the  calculated  vibrational  relaxation  rates.  Section  V  contains 
conclusions  and  a  discussion  of  several  issues  for  which  further 
work  is  needed  to  resolve. 
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II.  COMPUTATIONAL  METHODS 
A.  Potential  energy  surface 

The  nitromethane  intramolecular  potential  was  described 
by  the  Sorescu,  Rice,  and  Thompson2'27  (SRT)  force  field, 
which  consists  of  a  sum  of  Morse  functions  for  the  bond 
stretches,  harmonic  oscillators  for  bond  angles,  and  cosine 
series  for  the  HCNO  dihedral  angles  and  the  improper  dihedral 
angle  for  the  C-NO2  out-of-plane  deformation.  Nitromethane  - 
Ar  and  Ar-Ar  interactions  were  described  by  the  Buckingham 
potential 

Vali  =  Aafse~  ,  (1) 

where  R  is  the  distance  between  atoms  a  and  fi.  The  values 
of  the  Ar-Ar  parameters,  AAr-Ar  =  695,963.0  kcal/mol,  B,\r  ,\r 
=  0.24748  A,  and  Cat-ac  =  1294.81  A6  kcal/mol,  were  ob¬ 
tained  by  fitting  the  Colbourn  and  Douglas28  potential  to 
Eq.  (1).  The  Colbourn  and  Douglas  potential,  which  is  given 
at  discrete  R  values,  is  based  on  spectroscopic  data  for  the 
five  lowest  vibrational  states  of  Ar2,  with  corrections  for  long- 
range  forces  based  on  fitting  to  the  experimental  second  virial 
coefficient.  Relative  one-standard  deviation  uncertainties  in 
the  fitted  parameters  are  +10%,  +1%,  and  ±1%  for  A  Ar- at, 
R Ar-Ar.  and  Cat-At.  respectively.  The  parameters  for  the  nitro¬ 
methane- Ar  interactions  were  obtained  using  the  combination 
rules 

Affy8  =  AaaApp, 

^cr/3  {j^aa  +  Bpp)/2, 

Ca/3  =  CaaCpp,  (2) 

where  a  is  Ar  and  f3  is  C,  N,  O,  or  H.  The  values  of  the  homo¬ 
atom  parameters  App,  Bpp,  and  Cpp  for  the  C-C,  N-N,  O-O, 
and  H-H  interactions  were  taken  from  Sorescu  et  al.29 

The  SRT  nitromethane  force  field  produces  harmonic 
frequencies  in  good  accord  with  ab  initio  calculations.2'1  These 
frequencies  yield  an  harmonic  zero  point  energy  of  30.3 
kcal/mol,  comparable  to  the  classical  internal  energy  of  50 
kcal/mol  at  the  start  of  the  trajectories  studied  here.  However, 
this  zero  point  energy  is  not  available  for  transfer  to  the  bath 
gas  and  for  this  reason  is  excluded  from  the  trajectory  calcula¬ 
tions.  The  harmonic  frequencies  also  permit  an  estimation  of 
the  quantum  dynamical  vibrational  density  of  states  using  the 
Whitten  Rabinovitch  approximation.10  By  this  estimate,  the 
number  of  harmonic  vibrational  states  per  wavenumber  varied 
from  ~500  to  300,000  over  the  internal  vibrational  energy 
range  of  20.0  kcal/mole  to  45.0  kcal/mol  typically  sampled  by 
our  classical  trajectories.  With  such  a  large  density  of  states, 
quantum  effects  should  not  be  prominent  over  this  range  of 
energies.  At  much  lower  vibrational  energies,  quantum  effects 
can  become  more  important.  As  an  example,  at  300  K,  the  15 
vibrational  degrees  of  freedom  of  nitromethane  should  lead 
to  a  thermal  vibrational  energy  of  8.94  kcal/mol,  that  is,  15 
kT  where  k  is  the  Boltzmann  constant.  However,  the  quantum 
harmonic  energy  above  the  zero  point  energy  for  nitromethane 
is  only  1.00  kcal/mol.  '1  In  effect,  the  vibrational  quantum  of 
most  of  the  degrees  of  freedom  is  much  larger  than  kT,  leading 
to  a  quantized  form  of  nitromethane  that  behaves  as  if  it  is 


a  classical  molecule  with  much  fewer  degrees  of  freedom. 
Our  classical  trajectory  study  follows  in  the  tradition  of  many 
published  trajectory  studies  of  relaxation,  avoids  the  tail  end 
of  vibrational  relaxation  where  quantum  effects  can  become 
more  important,  and  should  provide  a  useful  description  of 
relaxation  processes  in  nitromethane. 

B.  Simulations 

Initial  conditions  for  the  relaxation  simulations  were  ob¬ 
tained  in  four  steps:  (1)  determination  of  the  equilibrium  vol¬ 
umes  at  300  K  for  the  eight  pressures  studied,  (2)  selection  of 
microcanonical  phase  space  points  for  the  internal  degrees  of 
freedom  of  the  excited  molecule,  (3)  equilibration  of  the  Ar 
bath  around  the  molecule,  with  the  molecule  held  rigid  and 
stationary  during  equilibration,  and  (4)  assignment  of  thermal 
translational  velocities  to  the  molecule. 

All  simulations  with  the  bath  gas  were  performed  using 
the  LAMMPS 32  code,  with  three-dimensional  periodic  bound¬ 
ary  conditions  for  a  simulation  box  with  initial  dimensions  of 
60.0  Ax 60.0  Ax 60.0  A  containing  one  molecule  and  1000 
Ar  atoms.  The  starting  geometry  of  the  molecule  corresponds 
to  the  SRT  equilibrium  geometry  and  the  center  of  mass  of  the 
molecule  was  initially  placed  at  the  center  of  the  box,  that  is, 
at  X  =  Y  —  Z  =  30.0  A.  The  initial  positions  of  the  Ar  atoms 
were  chosen  randomly  as 

X  =  fxL, 

Y  =  {«L, 

Z  =  {zL,  (3) 

where  //  (i  =  X,  Y,  Z )  is  a  pseudorandom  number  uniformly 
distributed  on  (0,1)  and  L  is  the  initial  edge  length  of  the  cubic 
simulation  cell.  If  the  distance  of  the  kth  Ar  atom  from  an 
atom  of  the  nitromethane  molecule  or  any  of  the  k  —  1  Ar 
atoms  already  in  the  simulation  cell  was  less  than  1 .0  A,  then  a 
new  position  was  chosen  (the  Buckingham  potential  of  Eq.  (1) 
passes  through  a  maximum  at  unrealistically  small  values  of 
R  [ranging  from  0.69  A  to  0.77  A  for  the  Ar-Ar,  Ar-C,  Ar¬ 
bi,  Ar-O,  and  Ar-H  pairs  used  here]  then  diverges  to  negative 
infinity,  necessitating  the  1.0  A  exclusion  test).  Starting  from 
these  system  configurations  200  ps  isobaric-isothermal  ( NPT ) 
simulations  were  run  at  T  =  300KandP=  10,50,75, 100, 125, 
150,  300,  and  400  atm  using  a  0.10  fs  time  step,  barostat  and 
thermostat  coupling  parameters  10.0  fs  and  100.0  fs,  respec¬ 
tively,  and  a  simulation  box  constrained  to  remain  cubic.  Initial 
velocities  for  all  Ar  and  nitromethane  atoms  were  selected 
from  the  Maxwell  distribution.  The  last  100  ps  of  the  NPT 
simulation  were  used  to  calculate  the  average  length  of  the 
simulation  cell  edge  used  for  each  pressure,  ( L )  =  159.75781 
A  for  10  atm,  92.84972  A  for  50  atm,  80.86864  A  for  75  atm, 
73.24327  A  for  100  atm,  67.96773  A  for  125  atm,  63.87542 
A  for  150  atm,  51.41449  A  for  300  atm,  and  47.77249  A  for 
400  atm.  These  edge  lengths  were  used  to  define  the  volumes 
of  the  cells  for  all  subsequent  isochoric-isothermal  ( NVT )  and 
isochoric-isoergic  (NVE)  simulations. 

Microcanonical  phase  space  points  for  the  isolated  excited 
molecule  were  obtained  using  the  GENDYN33  code.  Starting 
with  the  molecule  at  the  SRT  equilibrium  geometry,  a  Markov 
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walk34  of  two-million  steps  was  carried  out  in  configuration 
space  to  move  the  system  away  from  equilibrium  to  a  random 
geometry  consistent  with  the  microcanonical  ensemble  for 
the  internal  degrees  of  freedom  with  no  restriction  on  the 
angular  momentum.  At  each  Markov  step,  each  of  the  21 
Cartesian  coordinates  was  displaced  by  random  fractions  of 
+0.055  A  subject  to  the  restriction  that  none  of  the  six  covalent 
bond  lengths  in  the  molecule  became  larger  than  2.0  A;  these 
parameters  gave  an  acceptance/rejection  ratio  of  0.48  for  the 
molecular  internal  energy  studied  (50.0  kcal/mol).  The  initial 
phase  space  point  for  the  internal  degrees  of  freedom  of  the 
excited  molecule  for  the  first NVE  relaxation  trajectory  was  ob¬ 
tained  by  combining  the  configuration  at  the  end  of  the  initial 
random  walk  with  a  set  of  2 1  Cartesian  momenta  selected  from 
the  mass-weighted  normal  distribution  with  zero  mean  and 
unit  variance,  removing  the  center-of-mass  momentum,  and 
scaling  the  atomic  momenta  to  give  the  desired  microcanonical 
energy.  Molecular  initial  conditions  for  subsequent  trajectories 
were  similarly  obtained  based  on  configurations  sampled  at 
regular  10  000  step  intervals  along  the  Markov  chain. 

Following  the  selection  for  the  initial  conditions  for  the 
molecule,  the  molecular  center  of  mass  was  located  at  the 
center  of  the  simulation  box  with  edge  length  ( L(P )).  Then, 
the  molecule  was  held  fixed  and  rigid  as  new  positions  of  all 
Ar  atoms  were  chosen  randomly  using  Eq.  (3)  with  L  =  ( L(P )) 
and  new  velocities  for  the  Ar  atoms  were  selected  from  the 
Maxwell  distribution.  Nose-Hoover33  NVT  simulations  of  100 
ps  duration  were  performed  at  T  —  300  K  with  time  step  0.10 
fs  and  thermostat  coupling  parameter  10.0  fs.  After  the  ther- 
malization  of  the  bath,  the  internal  velocities  of  the  molecule 
were  restored  and  thermal  velocities  for  the  molecular  center 
of  mass  were  selected  from  the  Maxwell  distribution  at  300 
K.  This  completes  the  description  of  the  initial  conditions 
selection. 

Trajectories  were  calculated  in  the  NVE  ensemble,  with 
fixed  time  step  0.05  fs,  to  determine  the  relaxation  rates.  Tests 
were  performed  for  P  =  100  atm  for  systems  containing  500, 
1000,  and  1500  Ar  atoms.  It  was  found  that  the  decay  of  the 
excess  vibrational  energy  of  the  molecule  was  converged  with 
1000  Ar  atoms.  Ensembles  of  1000  trajectories  were  computed 
for  each  pressure.  Each  trajectory  was  integrated  for  1000  ps, 
except  for  10  atm  for  which  the  trajectories  were  integrated  for 
5000  ps. 

Translational,  vibrational,  and  rotational  energies  of  the 
nitromethane  molecule  were  calculated  as  functions  of  time. 
The  separation  between  vibrational  and  rotational  energies  of 
the  molecule  at  a  given  time  was  approximated  as  follows: 
After  removing  the  instantaneous  center-of-mass  velocity  and 
position,  the  internal  energy  of  the  molecule  was  defined  as 
the  sum  of  the  intramolecular  potential  energy  and  the  re¬ 
maining  kinetic  energy.  The  vibrational  energy  was  debited  as 
the  energy  of  the  molecule  that  remained  after  subtracting  the 
molecular  angular  momentum,  and  the  rotational  energy  was 
debned  as  the  difference  between  the  internal  energy  and  the 
vibrational  energy. 

As  mentioned  above,  at  pressures  greater  than  48  atm, 
Ar  is  a  supercritical  buid.  To  our  knowledge,  MD  calcula¬ 
tions  have  not  been  performed  with  the  potential  used  here  to 
model  Ar  in  the  supercritical  regime.  However,  a  quite  similar 


pairwise  additive  cib  initio  Ar-Ar  potential36  has  been  used  in 
MD  calculations  '7  of  500  Ar  atoms  with  periodic  boundary 
conditions  at  a  temperature  of  350  K  and  a  pressure  range  of 
~50  atm  to  ~900  atm.  The  calculated  atomic  pair  correlation 
function  agreed  well  with  the  experimental  one  measured  by 
neutron  diffraction.  ’4  This  suggests  that  the  simulations  re¬ 
ported  here  appropriately  describe  the  supercritical  regime. 

III.  RESULTS  AND  DISCUSSION 

We  have  calculated  trajectories  to  determine  the  relaxation 
of  nitromethane,  with  initial  internal  energy  50.0  kcal/mol 
above  the  classical  ground  state,  in  Ar  at  300  K  and  at  pressures 
of  10,  50,  75,  100,  125,  150,  300,  and  400  atm.  In  the  brst  part 
of  this  section,  the  ensemble-averaged  energy  in  nitromethane 
as  a  function  of  time  is  presented  for  all  eight  pressures,  and 
the  qualitative  features  of  the  decay  rates  are  discussed.  In  the 
second  part  of  this  section,  analytic  bts  to  the  decay  curves  are 
presented,  which  allows  for  a  more  precise  and  quantitative 
discussion  of  these  features. 

A.  Decay  curves 

The  ensemble-averaged  rotational  (blue  curve),  vibra¬ 
tional  (red  curve),  and  translational  (green  curve)  energies 
over  1000  ps  are  shown  in  Fig.  1  for  100  atm.  The  results 
clearly  show  the  equilibration  of  the  rotational  and  transla¬ 
tional  energy  to  values  approaching  1.5  kT  ( — 0.9  kcal/mol 
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FIG.  1.  Ensemble-averaged  vibrational,  rotational,  and  translational  energies 
versus  time  at  100  atm.  The  inset  is  a  zoom-in  of  the  lower  left-hand  corner 
of  the  plot. 
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for  T  =  300  K).  The  inset  in  Fig.  1  shows  that  rotation  and 
translation  are  equilibrated  within  ~15  ps,  while  the  full  plot 
shows  that  the  vibrational  energy  does  not  reach  equilibrium 
within  1000  ps.  These  equilibration  times  vary  with  pressure 
in  expected  ways  but  the  two-orders-of-magnitude  difference 
in  equilibration  times  for  rotation  and  vibration  stays  the 
same.  As  the  green  curve  in  the  inset  in  Fig.  1  shows,  the 
translational  energy  of  the  molecule  deviates  slightly  from 
the  equilibrium  value  for  the  first  few  picoseconds.  This  is  an 
artifact  of  the  initial  conditions  selection:  Interactions  between 
the  molecule  and  the  bath  were  neglected  during  assignment 
of  the  molecular  center-of-mass  kinetic  energy  whereas,  in 
reality,  the  forces  between  the  molecule  and  the  gas  atoms  in 
the  starting  configuration  will  not  generally  be  precisely  zero. 

In  Fig.  2,  the  normalized  ensemble-averaged  excess  vibra¬ 
tional  energy  Enrm.vu,{t)  is  plotted  as  a  function  of  time  for 
all  eight  pressures  in  the  main  panel  for  t  <  1000  ps  and  for 
P  —  10  atm  in  the  inset  for  t  <  5000  ps.  Excess  energy  at 
a  given  time  is  the  energy  above  the  thermal  energy  at  that 
time  including  the  temperature  increase  of  the  bath  gas  due  to 
molecular  energy  lost  up  to  that  time  (with  full  thermalization 
of  the  initial  molecular  excitation  of  50  kcal/mol  over  1000 
Ar  atoms,  the  temperature  increase  is  ~13  K).  The  results 
are  normalized  to  the  excess  energy  appropriate  to  the  initial 
absolute  value  of  the  energy  shown  in  Fig.  1,  which  is  the 
same  in  all  cases.  The  excess  vibrational  energy  is  plotted  on  a 
log  scale,  for  which  single-exponential  decay  would  be  linear 
(because  the  asymptotic  value  of  the  excess  energy  is  zero). 
The  results  show  curvature  at  the  higher  pressures  indicating 


FIG.  2.  Enrm.vib  versus  time.  The  solid  curves  are  the  trajectory  results  and 
the  dashed  curves  are  fits  to  the  trajectory  results.  The  inset  is  an  expansion 
of  the  10  atm  results  to  5000  ps. 


multiple  decay  rates  over  the  temporal  interval  of  the  plot.  For 
the  lowest  pressure,  10  atm,  it  appears  that  a  single  decay  rate 
applies  over  1000  ps;  however,  as  the  results  in  the  inset  show, 
when  the  process  is  followed  for  5000  ps,  there  is  clearly  a 
curvature.  The  results  in  Fig.  2  show  that  there  is  not  a  single 
decay  rate  at  any  pressure  as  the  vibrational  energy  relaxes 
to  equilibrium.  The  results  in  Fig.  2  also  clearly  show  that 
the  excess  vibrational  energy  decays  at  increasing  rates  with 
increasing  pressure. 

It  is  useful  to  examine  the  behavior  of  the  rotational  relax¬ 
ation  on  a  much  shorter  time  scale  than  shown  in  Fig.  2 
for  vibrational  relaxation.  The  results  shown  in  Fig.  3  are 
for  the  excess  rotational  energy  Enrm.I0t(t)  for  the  first  50  ps 
except  for  10  atm,  which  is  shown  for  200  ps  (see  inset).  The 
most  important  feature  of  Fig.  3  is  that,  as  was  the  case  for 
the  vibrational  relaxation  shown  in  Fig.  2,  excess  rotational 
energy  decays  with  a  rate  that  increases  monotonically  with 
pressure.  However,  the  most  striking  feature  of  the  results  in 
Fig.  3  is  that  approximately  5%  of  the  initial  excess  rotational 
energy  remains  in  the  molecule.  Although  not  shown  in  the 
figure,  at  1000  ps,  some  of  this  5%  of  the  initial  rotational 
energy  still  remains.  As  seen  in  the  inset  of  Fig.  1,  5%  of 
the  initial  rotational  energy  amounts  to  a  couple  of  tenths 
of  a  kcal/mol.  Some  of  this  residual  excess  energy  could  be 
due  to  the  approximate  nature  of  separation  of  rotational  and 
vibrational  energies.  Also,  the  statistical  variations  due  to  finite 
sample  size  are  comparable  to  this  slow  decay  feature  as  seen 


FIG.  3.  Enrm. rot  versus  time  analogous  to  Fig.  2  only  for  a  maximum  time 
of  50  ps  (200  ps  for  the  inset).  Black  shading  in  the  inset  is  the  statistical 
uncertainty. 
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in  the  black  shading  in  the  inset,  which  represents  one-cr 
statistical  error  in  the  mean  for  10  atm  (cr  values  at  other 
pressures  are  similar).  Finite  sampling  statistics  also  affect  the 
value  of  excess  rotational  energy  through  statistical  variations 
in  the  rising  temperature  of  the  bath  gas.  To  understand  more 
about  this  slowly  decaying  residual  of  rotational  energy,  the 
ensemble  size  would  need  to  be  significantly  larger,  an  effort 
that  has  not  been  undertaken  for  this  paper. 

B.  Analytic  fits  to  decay  curves 

The  presence  of  curvature  in  Fig.  2  for  the  excess  vibra¬ 
tional  energy  requires  a  representation  more  flexible  than 
a  single  exponential.  For  example,  Paul  et  a/.23  used  a  bi¬ 
exponential  fit  to  represent  curvature  in  their  QU^/No  results. 
While  this  form  can  fit  our  data  reasonably  well,  it  has  two 
disadvantages.  First,  if  one  requires  the  two  coefficients,  one 
for  each  exponential,  to  sum  to  the  initial  value  of  the  energy, 
there  are  three  adjustable  parameters.  For  some  of  the  data  in 
Fig.  2,  there  are  strong  correlations  between  the  three  param¬ 
eters,  leading  to  large  uncertainties  in  the  three  optimized 
values.  Second,  the  two  decay  rates  in  some  cases  are  not 
sufficiently  different  for  one  rate  to  completely  dominate  in 
some  portion  of  the  temporal  interval  of  Figs.  2  and  3,  which 
complicates  separate  interpretations  of  each  rate.  Another  way 
to  say  this  is  that  there  is  no  physical  motivation  for  just  two 
rates.  Rather,  the  decay  rate  is  continuously  changing  as  the 
energy  of  the  molecule  transfers  to  the  bath  gas  and  the  fits  do 
not  capture  this  with  just  two  rates. 

An  alternative  function  was  introduced  by  Lendvay  and 
Schatz19  (LS)  who  integrated  out  the  often  used  relationship 
A E  oc  E(t)m  from  isolated-molecule/bath-gas  collision  studies, 
where  A E  is  the  average  energy  lost  per  collision  by  the  mole¬ 
cule  with  initial  energy  E(t).  Because  the  A E  is  proportional 
to  the  derivative  of  E(t ),  A E  oc  E(t)m  can  be  converted  into  a 
differential  equation  whose  integration  gives 

^r=Enrm(t)  =  [l-(l-m)kit]i^.  (4) 

E(  0) 

This  form  for  the  decay  of  the  normalized  energy  Emm(t)  has 
the  virtue  that  it  involves  only  two  parameters,  not  three,  and 
each  parameter  more  clearly  responds  to  different  temporal 
regions  of  the  decay.  At  t  =  0,  the  k\  parameter  equals  dEnn„/dt, 
the  initial  rate  of  decay.  To  the  extent  that  this  form  is  appro¬ 
priate,  k\  will  be  determined  primarily  by  the  early  behavior 


of  the  decay  curve.  The  parameter  m  determines  curvature. 
When  m  =  1,  Eq.  (4)  can  be  shown  '8  to  collapse  to  a  single 
exponential  function  with  rate  kv  For  m  >  1 ,  Emm(t)  curves  up 
from  the  single-exponential  straight-line  behavior  on  a  semi¬ 
log  plot,  meaning  that  as  energy  flows  out  of  the  molecule,  the 
rate  of  loss  of  the  remaining  energy  slows  down.  For  m  <  1 
(including  negative  values),  Enn„(t )  curves  down  from  single¬ 
exponential  (linear)  behavior,  meaning  that  as  energy  flows  out 
of  the  molecule,  the  rate  of  loss  of  the  remaining  energy  accel¬ 
erates.  While  m  <  1  can  result  in  an  optimal  fit  of  Eq.  (4)  to 
ensemble-averaged  energy  decay  over  a  finite  interval  of  time, 
it  is  not  physically  meaningful  because  the  term  inside  the 
brackets  in  Eq.  (4)  will  become  negative  at  large  values  of  time, 
producing  an  Enrm(t )  that  is  no  longer  a  real  number.  Because 
Eq.  (4)  cannot  be  derived  from  more  rigorous  collision  theory, 
the  expression  is  only  a  fitting  function,  however,  one  with  only 
two  parameters,  each  of  which  has  a  useful  phenomenological 
interpretation. 

The  application  of  the  LS  fit  for  excess  vibrational  energy 
decay  is  represented  as  dashed  curves  in  Fig.  2.  The  optimal 
parameter  values  for  m  and  k\  are  listed  in  Table  I  along  with 
the  rms  absolute  error  for  the  normalized  decay  £’„rai-vih(0- 
Also  listed  are  the  uncertainties  Am  and  Ak,  determined  by  a 
bootstrap  method  9  using  2000  random  realizations  of  the  1000 
trajectories  calculated  at  a  given  pressure.  For  all  pressures,  the 
fit  has  an  rms  error  less  than  0.01.  Inspection  of  Fig.  2  shows 
that  while  the  fit  is  quite  good  for  each  pressure,  at  higher 
pressures  the  computed  Enrm.v^(t)  becomes  more  ragged  at 
longer  times.  This  raggedness  is  probably  due  to  the  finite 
sampling  error  becoming  comparable  to  the  quite  small  excess 
vibrational  energies  at  these  long  times.  Although  it  is  not 
apparent  at  the  scale  of  Fig.  2,  the  LS  fit  for  the  10  atm  results 
is  somewhat  inaccurate  at  early  times;  this  is  not  typical  of  the 
fits  at  other  nearby  pressures.  Figure  4  is  a  zoom-in  of  the  first 
400  ps  of  Fig.  2  for  the  five  lowest  pressures.  Unlike  the  results 
at  all  the  other  pressures,  Enrm.v]b(t)  at  10  atm  initially  curves 
downward.  This  would  require  a  value  of  m  less  than  1  but 
the  inset  in  Fig.  2  shows  clearly  that  over  5000  ps  Enrm.^(t) 
does  substantially  curve  upward,  leading  to  least-squares  value 
m  =  1.34  (see  Table  I).  Figure  5  shows  plots  of  the  values  of 
m  and  k,  as  functions  of  the  time  interval  of  the  10  atm  data 
included  in  the  fit,  always  starting  from  t  -  0.  Clearly,  LS  fits 
confined  to  early  times  have  values  of  m  <  1  and  k\  values 
~80%  of  the  values  listed  in  Table  I.  The  implication  of  this  is 
that  the  least-squares  value  of  k,  for  10  atm  (see  Table  I)  should 


TABLE  I.  As  functions  of  pressure  and  density,  parameter  values  m  and  k\  for  vibrational  relaxation,  their 
uncertainties  Am  and  A/:j,  and  the  rms  fitting  error  for  the  LS  fit  to  E^r/w-vibW- 


Pressure 

(atm) 

Density  (xlO3  A  3) 

m 

Am 

h,  (xlO3  ps^1) 

A/q(xl0sps  *) 

rms  error  (xlO3) 

10 

0.245 

1.34 

0.004 

0.341 

0.09 

3.11 

50 

1.25 

1.46 

0.045 

1.04 

1.92 

1.25 

75 

1.89 

1.58 

0.037 

1.43 

2.58 

1.89 

100 

2.55 

1.56 

0.031 

1.75 

2.96 

2.55 

125 

3.19 

1.53 

0.027 

2.00 

3.09 

3.18 

150 

3.84 

1.60 

0.024 

2.41 

3.69 

3.84 

300 

7.36 

1.42 

0.018 

3.68 

5.38 

7.36 

400 

9.17 

1.37 

0.018 

4.40 

6.54 

9.17 
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FIG.  4.  As  in  Fig.  2,  except  a  zoom-in  of  the  first  400  ps  and  only  for  the  five 
lowest  pressures. 

be  interpreted  as  a  more  global  value  of  the  rate  rather  than  the 
actual  initial  relaxation  rate,  which  (see  Fig.  5)  is  noticeably 
slower.  As  Fig.  4  shows,  this  difficulty  is  unique  to  the  10  atm 
results  and  may  represent  some  small-scale  imperfections  in 


end  time  (ps) 


* 

o 

o 

o 

o 


FIG.  5.  Optimized  values  of  m  and  It,  for  vibrational  relaxation  at  10  atm 
plotted  versus  the  time  interval  of  the  Enrm-vib  data.  Parameters  are  color  coded 
to  the  appropriate  axis.  Dashed  lines  are  straight  line  segments  to  aid  the  eye. 


the  preparation  of  the  initial  conditions  of  the  simulation.  This 
warrants  further  study. 

At  lower  pressures  corresponding  more  closely  to  ideal 
behavior,  the  density  and  pressure  are  directly  proportional 
to  each  other;  however,  at  higher  pressures,  the  density  more 
slowly  increases  in  response  to  pressure.  The  first  noticeable 
break  in  direct  proportionality  can  be  seen  in  the  last  few  rows 
of  Table  I.  Extreme  pressures  are  necessary  to  condense  the  Ar 
gas  to  the  solid  state  at  300  K,  but  it  is  instructive  to  investigate 
how  much  the  density  for  the  highest  value  studied  here  (see 
Table  I)  differs  from  the  nominal  solid-state  density.  Using 
the  Ar-Ar  potential  described  in  Sec.  II  A,  the  Ar  hard-sphere 
radius  is  1.675  A.  Given  that  radius  as  the  spherical  size  of  Ar, 
the  packed-sphere  formula40  of  Gauss  provides  the  maximum 
number  of  spheres  that  can  be  packed  into  a  given  volume. 
That  density  is  0.0374  A-3.  This  value  can  be  compared  to 
the  measured41  density  of  0.0275  A-3  for  face-centered-cubic 
crystal  Ar  at  4.2  K  and  saturated  vapor  pressure.  The  crystal 
structure  expands  with  increasing  temperature  and  contracts 
with  increasing  pressure42  for  an  unknown  end  result  at  300  K 
and  many  thousands  of  atmospheres  of  pressure.  Most  likely, 
the  packed-sphere  density  0.0374  A-3  is  an  upper  bound  for  the 
true  density  of  the  solid  phase  of  Ar  at  300  K  and  by  that  guide 
the  maximum  density  studied  here  is  ~25%  of  the  solid-state 
density. 

Plots  of  the  parameter  values  m  and  k,  as  functions  of 
density  are  shown,  respectively,  in  Figs.  6  and  7.  The  results  in 
Fig.  6  show  that  the  m  parameter  assumes  values  in  a  narrow 
range  of  1.3-1. 6  over  a  factor  of  ~40  variation  in  the  density. 


FIG.  6.  From  Table  1,  m  values  for  vibrational  relaxation  with  their  uncer¬ 
tainty  Am  versus  density.  The  solid  curve  is  a  fit  (see  text  for  details). 
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FIG.  7.  From  Table  I,  k,  values  for  vibrational  relaxation  with  their  uncer¬ 
tainty  Ah,  versus  density.  The  solid  curve  is  a  fit  (see  text  for  details).  The 
single  open  circle  is  the  early  time  value  of  k,.  taken  from  Fig.  5. 

The  persistence  of  this  band  of  curvature  over  such  a  wide 
variation  of  conditions  suggests  that  the  unchanging  internal 
structure  of  nitromethane  is  involved.  Perhaps  the  effective 
coupling  of  high-  to  low-frequency  vibrational  modes  changes 
in  the  molecule  as  energy  is  lost,  or  perhaps  the  degree  of 
coupling  of  vibrational  and  rotational  energies  changes  with 
energy  loss.  Such  changes  could  influence  how  readily  the 
molecular  internal  energy  is  lost  to  the  bath  and  result  in 
a  slowdown  of  the  rate  as  the  energy  decreases.  Additional 
studies  to  uncover  the  mechanisms  responsible  for  curvature 
are  planned.  The  CgF^Ar  study  of  Paul  et  al.23  also  exhibits 
curvature  (although  not  characterized  by  a  parameter  m)  at  all 
pressures  studied. 

The  solid  curve  in  Fig.  6  is  a  fit  of  the  density  dependence 
of  m  to  the  functional  form  a+bpe~cp ,  where  p  is  the  density. 
Because  the  mechanisms  of  curvature  are  not  yet  known,  this 
functional  form  is  one  of  convenience  and  is  not  based  on  any 
physical  justifications.  For  the  least-squares  parameter  values 
( a,b,c )  =  ( 1 .252,292.7A\338.8A3),  this  functional  form  repre¬ 
sents  the  results  within  a  two-sigma  uncertainty  with  Am  in 
Table  I  being  sigma.  The  least-squares  value  of  a  suggests  that 
substantial  curvature  exists  for  all  gas-phase  densities. 

The  results  in  Fig.  7  show  that  the  k,  parameter  assumes 
values  that  increase  monotonically  over  the  factor  of  -40 
variation  in  the  density.  However,  the  results  in  the  figure  also 
clearly  show  that  k\  is  not  directly  proportional  to  density  over 
this  interval.  At  the  density  of  -0.002  A'3,  corresponding  to 
-100  atm,  k,  begins  to  drop  below  what  is  expected  from  a 
direct  proportionality.  Recall  that  k,  is  the  relaxation  rate  at 


t  =  0  for  the  LS  fit.  The  initial  conditions  of  the  molecule 
are  statistically  identical  for  all  densities  (ignoring  the  small 
artifact  noted  for  the  translational  energy  in  the  discussion  of 
Fig.  1).  This  implies  that  the  loss  of  direct  proportionality  is  a 
consequence  of  a  change  in  the  nature  of  the  interaction  of  Ar 
bath  with  the  molecule  as  a  function  of  density.  A  proposed 
mechanism  for  this  change,  described  in  Sec.  IV,  leads  to  a 
functional  form  whose  least-squares  optimization  leads  to  the 
solid  curve  in  Fig.  7. 

The  LS  fits  for  excess  rotational  energy  decay  are  shown 
by  dashed  curves  in  Fig.  3.  However,  the  presence  of  the  slowly 
decaying  residual  energy  makes  this  fitting  not  straightfor¬ 
ward.  Fitting  rotational  decay  over  the  full  1000  ps  of  the 
simulation  sacrifices  an  accurate  representation  of  the  major 
early  decay  for  a  more  reliable  representation  of  the  minor 
residual  decay.  Given  that  the  mechanism  of  the  residual  decay 
is  not  known,  a  physically  motivated  functional  form  different 
from  and  more  flexible  than  LS  cannot  now  be  developed. 
Rather,  the  approach  taken  for  the  results  in  Fig.  3  is  to  limit 
at  each  pressure  the  temporal  interval  over  which  the  LS  fit 
is  applied.  The  error  of  the  mean  due  to  finite  sampling  is 
available  at  each  value  of  Enrm.m(t)  for  each  pressure  in  Fig. 
3  (see  the  black  shading  in  the  inset),  and  the  limit  selected  is 
the  maximum  time  for  which  the  LS  fit  falls  within  this  error 
85%  of  the  time.  Obviously  85%  is  not  physically  motivated.  A 
higher  percentage  results  in  a  fit  over  a  shorter  span  that  leaves 
out  some  of  the  major  decay.  A  lower  percentage  results  in  a 
compromised  fit  that  includes  the  residual  tail.  Inspection  of 
Fig.  3  reveals  that  the  final  LS  fit  represents  the  major  decay 
well  with  the  possible  exception  of  £nmi-rot(f )  at  75  atm,  which 
exhibits  a  more  gradual  approach  to  the  residual  tail  than  do 
the  results  at  nearby  pressures. 

Analogous  to  Table  I,  the  optimal  parameter  values  for  m 
and  k\  are  listed  in  Table  II  along  with  the  rms  absolute  error  for 
the  normalized  decay  Enrm.rol(t)  and  the  maximum  time  over 
which  the  fit  was  applied.  The  fitting  error  is  -2  to  -35  times 
larger  than  that  for  £„™-vib(0  and  for  a  temporal  interval  that 
is  0.005-5  times  smaller.  This  is  obviously  a  reflection  of  the 
noisier  data  in  Fig.  3  due  to  the  relatively  small  magnitude 
of  Ennn-mxU )  compared  to  Enrm-V,b(t).  The  maximum  time  over 
which  the  LS  fit  is  applied  increases  monotonically  with  pres¬ 
sure  except  for  75  atm  for  which  the  approach  of  the  major 
decay  to  the  residual  decay  is  noticeably  slower.  Because  of 
the  less-than-definitive  nature  of  the  fit,  no  attempt  was  made 
to  determine  the  uncertainties  Am  and  Ak,.  However,  it  is  quite 
clear  that  the  m  value  for  rotational  relaxation  (Table  II)  is 
consistent  with  the  1 .3—1 .6  range  for  vibrational  relaxation 
(Table  I).  The  k,  value  for  rotational  relaxation  (Table  II)  is  at 
least  a  hundred  times  larger  than  the  corresponding  value  for 
vibrational  relaxation  (Table  I). 

A  two-sided  plot  of  both  the  m  and  k,  parameter  values 
for  rotational  relaxation  from  Table  II  as  functions  of  density 
is  shown  in  Fig.  8.  As  with  the  vibrational  m  dependence,  the 
rotational  m  dependence  on  density  in  Fig.  8  is  relatively  small 
at  the  low  end  and  high  end  of  the  density  interval  studied  and 
exhibits  a  maximum  for  intermediate  density  values.  Unlike 
the  vibrational  k,  dependence  (Fig.  7),  the  k,  dependence  on 
density  in  Fig.  8  is  directly  proportional  to  density  as  indi¬ 
cated  by  the  least-squares  straight  line  with  a  zero  intercept. 
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TABLE  II.  As  a  function  of  pressure  and  density,  parameter  values  m  and  In  for  rotational  relaxation,  the  rms 
fitting  error,  and  the  maximum  time  for  the  LS  fit  to  Enrm  mU ). 


Pressure  (aim) 

Density  (xlO-3A~3) 

m 

kid*-1) 

rms  error  (x  10“  *) 

Max  time  (ps) 

10 

0.245 

1.42 

0.030 

0.87 

195.4 

50 

1.249 

1.40 

0.177 

1.23 

27.8 

75 

1.891 

1.37 

0.286 

1.28 

14.0 

100 

2.545 

1.59 

0.466 

1.04 

21.8 

125 

3.185 

1.58 

0.537 

0.92 

19.4 

150 

3.837 

1.52 

0.663 

1.13 

7.9 

300 

7.358 

1.56 

1.50 

0.88 

7.0 

400 

9.172 

1.43 

1.68 

1.31 

4.9 

This  qualitative  difference  between  vibrational  and  rotational 
dynamics  is  discussed  in  Sec.  [V. 

IV.  A  PROPOSED  MODEL 

In  Figs.  6  and  7,  Enrm.v jb(f)  has  been  characterized  in 
terms  of  k,  and  m.  The  m  parameter  represents  the  curva¬ 
ture  over  the  entire  temporal  interval  of  the  simulations  and 
undoubtedly  involves  the  internal  coupling  of  modes  of  motion 
in  nitromethane  as  the  internal  energy  decreases  as  well  as 
die  changing  character  of  nitromethane/Ar  interactions  for 
different  densities.  The  curve  in  Fig.  6  is  the  density  depen¬ 
dence  fit  to  a  convenient  functional  form  unmotivated  by  any 
insight  into  the  molecular  mechanisms  responsible  for  curva¬ 
ture.  In  contrast,  the  k,  parameter  is  nominally  the  initial  rate 
of  decay  of  E„m  v,b(r).  Because  the  initial  conditions  of  the 
molecule  are  statistically  the  same  at  every  pressure,  the  pres¬ 
sure  dependence  of  k,  is  primarily  dependent  on  the  changes 


density  (1/A3 ) 


FIG  8.  From  Table  II.  m  and  k i  values  for  rotational  relaxation  versus 
density.  The  dashed  line  segments  connect  m  values.  The  solid  line  is  a  direct 
proportionality  fit  to  k\. 


in  the  nitromethane/Ar  interactions  as  a  function  of  density.  In 
this  section,  a  physical  model  will  be  developed  for  the  density 
dependence  of  the  nitromethane/Ar  interactions  that  motivates 
the  functional  form  used  to  obtain  the  fit  to  the  k,  pressure 
dependence  shown  in  Fig.  7. 

At  low  densities  for  which  the  isolated  binary  collision 
assumption  is  valid,  k,  should  be  directly  proportional  to 
density.  As  the  density  increases,  one  obvious  possible  change 
is  the  development  of  Ar  clusters.  If  the  clusters  are  numerous 
enough,  and  if  the  vibrational  relaxation  efficacy  of  a  nitro- 
methane/Arn  collision  is  sufficiently  different  from  that  of 
a  nitromethane/ Ar-atom  collision,  the  loss  of  linear  propor¬ 
tionality  between  k,  and  density  seen  in  Fig.  7  could  be  ex¬ 
plained.  A  simple  test  of  this  hypothesis  is  to  run  the  simulation 
without  any  van  der  Waals  wells.  To  this  end.  simulations  at 
10  atm  and  150  atm4'  were  performed  in  which  all  Ar-Ar 
and  nitromethane-Ar  pairwise  Buckingham  potentials,  Eq.  ( 1 ), 
were  truncated  at  the  bottom  of  each  well  and  shifted  vertically 
such  that  the  energy  at  and  beyond  the  truncation  point  is 
zero.  This  means  that  all  interactions  of  Ar  with  itself  or  nitro¬ 
methane  are  purely  repulsive  and  therefore  no  Ar  clusters  or  Ar- 
molecule  complexes  can  be  formed.  The  resulting  fwvibfr) 
was  fit  to  the  LS  function  and  yielded  k\  values  that  are  49%  (for 
10  atm)  and  47%  (for  150  atm)  as  large  as  the  corresponding 
values  shown  in  Table  I.  The  fact  that  these  two  percentages 
are  so  close  in  value  indicates  that  the  scale  of  the  functional 
dependence  on  pressure  has  changed  but  the  shape  has  not. 
This  means  that  in  the  absence  of  van  der  Waals  clusters  of  any 
sort,  k j  still  departs  from  proportionality  to  density  at  higher 
densities.  The  fact  that  kj  is  smaller  in  the  absence  of  potential 
wells  is  generally  consistent  with  the  physical  notion44  that 
wells  increase  the  relative  velocity  of  the  collision  just  prior 
to  contact  with  the  repulsive  wall  of  the  molecule. 

The  cluster  hypothesis  can  be  regarded  as  a  restricted  form 
of  a  more  general  model  wherein  increasing  density  results 
in  the  molecule  interacting  simultaneously  with  multiple  Ar 
atoms.  At  condensed-phase  densities,  the  molecule  would  be  in 
contact  with  Ar  atoms  over  its  entire  surface  area.  As  discussed 
in  Sec.  III.  at  400  atm  (the  highest  pressure  studied)  the  density 
is  about  25%  of  the  nominal  solid-state  density.  Of  course,  at 
low  enough  densities,  nitromethane  will  only  interact  with  one, 
if  any,  Ar  atom  at  a  given  time.  As  the  density  is  increased,  the 
probability  of  a  nitromethane  collision  event  involving  two  or 
more  Ar  atoms  simultaneously  increases.  If  the  per-collision 
efficacy  of  Enm-vib(f)  relaxation  varies  with  the  number  of  Ar 
atoms  in  a  single  collision  event,  then  the  loss  of  direct  propor- 
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tionality  of  seen  in  Fig.  7  could  be  explained.  Note  that  this 
explanation  depends  on  multiple,  simultaneous  molecule/bath- 
gas  interactions  due  to  high-density  proximity.  As  such  it  is 
consistent  with  bound  interactions  as  in  clusters  and  complexes 
but  does  not  require  that  the  interactions  be  bound. 

If  one  approximates  both  Ar  and  nitromethane  as  hard 
spheres,  one  can  render  the  notional  model  outlined  above 
into  an  exercise  in  combinatorial  arithmetic.  Imagine  plac¬ 
ing  the  nitromethane  hard  sphere  at  the  center  of  a  large  but 
finite  volume  and  discretizing  all  of  the  volume  outside  of 
nitromethane  into  cells  that  can  hold  at  most  one  Ar  hard 
sphere.  Figure  9  is  a  two-dimensional  illustration  of  this  model. 
There  will  be  a  finite  number  of  cells  that  are  adjacent  to  the 
nitromethane  (as  indicated  by  the  white  squares  in  Fig.  9) 
but  the  vast  majority  of  the  cells  do  not  border  nitromethane 
(yellow  squares  in  Fig.  9).  A  given  simulation  density  sets  how 
many  Ar  atoms  must  go  into  the  finite  volume  to  be  consistent 
with  that  density.  Because  of  the  hard-sphere  approximation, 
any  cell  in  the  volume  is  equally  likely  to  contain  an  Ar  atom 
but  no  cell  can  hold  more  than  one  Ar  atom.  Combinato¬ 
rial  arithmetic  can  determine  how  many  total  combinations 
of  Ar  placements  within  the  cells  there  are;  and  how  many 
combinations  have  no  Ar  atoms  touching  nitromethane  (Fig. 
9(a)),  one  Ar  atom  touching  nitromethane  (Fig.  9(b)),  two  Ar 
atoms  touching  nitromethane  (Fig.  9(c)),  etc.  This  leads  to 
an  analytic  density-dependent  distribution  of  discrete  collision 
events  involving  multiple  Ar  atoms  for  the  finite  volume  system 
from  which  the  limit  of  infinite  volume  can  be  obtained  in 
straightforward  fashion.  Within  the  context  of  the  hard-sphere 
approximation,  consider  a  snapshot  of  the  locations  of  the  Ar 
atoms  with  respect  to  nitromethane  at  the  same  given  time 
for  each  of  the  1000  trajectories  comprising  the  ensemble  at 
a  given  density.  For  the  given  time,  this  would  produce  an 
album  of  1000  pictures.  The  fraction  of  pictures  in  the  album 
where  one,  two,  three,  or  more  Ar  atoms  are  near  nitromethane 
should  be  approximated  by  this  combinatorial  model.  In  what 
follows,  the  mathematics  of  this  model  will  be  outlined  and 
the  parameters  necessary  to  apply  it  to  the  simulations  will 
be  discussed.  The  remainder  of  the  section  will  then  explore 
the  dynamical  implications  of  the  resulting  density-dependent 
distributions. 

Let  N  be  the  number  of  Ar-size  cells  outside  of  nitro¬ 
methane  in  a  large  but  finite  volume  that  encloses  the  nitro¬ 
methane.  Let  K  be  the  number  of  those  cells  that  “touch” 
nitromethane.  Let  M  be  the  number  of  Ar  atoms.  Let  Va,  be 
the  volume  of  one  Ar  atom,  that  is,  the  volume  of  one  cell.  If  p 
is  the  number  density  of  Ar  atoms  in  the  particular  simulation 
under  consideration,  then 

M  =  pVAlN.  (5) 

The  total  number  of  combinations  of  M  atoms  in  N  cells  is 

Nl 

total  combinations  = - .  (6) 

M\(N-  M) ! 


L  restricted  combinations 

K\  (, N-K)\ 

~  L\(K-L)\  (M  -  L)\(N  -  K  -  M  +  L)\' 


(7) 


The  total  number  of  combinations  of  L  atoms  in  K  cells  and 
M-L  atoms  in  N-K  cells  is 


FIG.  9.  2D  schematic  of  three  combinations  in  the  combinatorial  model  for 
molecule/multi- Ar  collisions.  White  squares  “touch”  red  CH3NO2. 
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In  the  limit  of  infinite  N 

..  (N  —K)\ 

Let  P(L)  be  the  normalized  probability  of  having  L  atoms  in 
K  cells,  that  is,  the  ratio  of  Eq.  (8)  to  Eq.  (6) 

p(L)=  U(^_Ly}PV^)L0-P^f-L  L  =  [0,1,..., AT],  (9) 

Inspection  of  this  result  shows  that  each  P(L)  matches  a  term  in 
the  binomial  expansion  of  (pV^+  1  -pVA[)K ,  which  of  course 
is  1,  confirming  a  normalized  probability. 

To  evaluate  P(L)  in  Eq.  (9),  VAr  and  K  must  be  determined. 
As  discussed  in  Sec.  Ill  B,  from  the  Ar-Ar  potential,  the  hard- 
sphere  radius  of  Ar  is  1 .675  A.  Vat  is  the  volume  of  the  sphere 
with  this  radius  and  K  is  the  maximum  number  of  Ar  atoms 
that  can  touch  the  surface  of  the  nitrome  thane  molecule,  which 
has  radius  /?ch3no2-  This  can  be  approximated  by  the  surface 
area  of  a  sphere  of  radius  Pch3no2  +  R\i  divided  by  the  cross- 
sectional  area  of  Ar.  Because  a  two  dimensional  array  of  circles 
that  touch  one  another  have  open  area  in  between  them,  the  Ar 
circular  cross-sectional  area  cannot  completely  cover  the  area 
of  the  sphere.  In  this  context,  a  more  appropriate  estimate  of 
the  cross-sectional  area  of  Ar  is  to  treat  Ar  as  a  square  with 
edge  length  2/?Ar.  The  result  is 


K  =  n^ 


1  + 


flcHjNOji2 

— ‘ 


(10) 


For  /?ch3no2.  a  value  of  2.2  A  is  appropriate  based  on  the 
equilibrium  geometry  of  nitromethane.15  Given  the  values  of 
R\t  and  /?ch3no2,  both  and  K  can  be  determined,  which  in 
turn  finalizes  for  nitromethane/Ar  the  P(L)  distribution  once 
the  density  p  is  specified. 

To  understand  the  dynamical  implications  of  this  distri¬ 
bution,  consider  the  simplified  relationship  between  a  rate 
constant  k  and  its  scattering  cross  section  cr  at  low  densities 


k  =  pcrv. 


(ID 


where  v  is  the  averaged  thermal  relative  velocity.  In  this  equa¬ 
tion,  cr  is  the  effective  cross  section  that  makes  crv  equiva¬ 
lent  to  the  integration  of  the  actual  collision-energy -dependent 
cross  section  over  the  Boltzmann  thermal  velocity  distribution. 
Equation  (11)  expresses  the  direct  proportionality  between  k, 
and  p  at  low  densities.  However,  this  equation  presumes  that 
there  is  only  one  dynamical  process,  namely,  the  collision  of 
one  bath-gas  atom  with  one  target  molecule. 

In  Eq.  (II),  pn  is  a  flux  of  atoms,  that  is,  it  is  the  number 
of  atoms  per  unit  time  per  unit  area  passing  through  an  arbi¬ 
trary  plane  in  space.  Using  P(L),  this  flux  can  be  decomposed 
into  fluxes  for  different  processes,  each  one  of  which  can,  in 
principle,  have  its  own  cross  section.  Because  P(L )  involves  L 
atoms,  the  density  can  be  decomposed  by 


=  pYj  yKP^P,j\  ~  3  Z  PL*  (12> 


tiU=iLP(D 


L=  1 


L=\ 


where  Pp(L)  is  normalized  to  unity  and  pl  is  the  number  of 
atoms  per  unit  volume  in  a  collision  event  involving  L  bath 


atoms  and  the  molecule  simultaneously.  This  decomposition 
leads 

K  K 

peer— »  YJPLVO-L  =  J^Pp(L)pvcrL.  (13) 

L= 1  L= 1 

Here,  crL  is  the  cross  section  per  atom  involved  in  the  collision 
event  and  therefore  the  overall  cross  section  for  the  process  is 
Lctl  .  This  follows  from  the  fact  that  pt  is  atom-based  while  the 
cross  section  is  for  a  flux  in  units  of  L  atoms.  Figure  10  shows 
Pp(L)  for  the  eight  densities  studied.  As  the  density  increases, 
the  distribution  shifts  to  multiple  Ar  collision  events. 

The  fit  to  k\  in  Fig.  7  is  a  least-squares  fit  of  ucrL  in 
Eq.  (13)  under  the  simplification  that  vctl  =  vcr2  for  all  L>  2. 
This  leaves  Eq.  (13)  with  only  two  adjustable  parameters:  txxi 
and  vcr2.  The  resulting  fit  is  quite  good,  with  an  rms  error  of 
5.3  x  10“5  ps_I  for  an  average  k,  value  of  2. 1  x  10-3  ps_! .  While, 
in  principle,  more  adjustable  parameters  can  be  used,  adding 
more  of  them  will  not  materially  improve  the  fit  and  indeed, 
for  a  sufficiently  large  number,  some  of  them  will  optimize  to 
aphysical  negative  values.  The  two-parameter  fit  has  a  qualita¬ 
tive  feature  consistent  with  the  k,  reported  in  Table  I  at  higher 
densities,  as  can  be  seen  from  the  following  rearrangement  of 
Eq.  (13): 

K 

pver  -»  Pp(\)pv{o-]-a-2)+YuPP(L)PVO'i 

L=  1 

=  Pp(l)pv(o-i-o-2)+ pvo-2.  (14) 

Specifically,  as  Pp(l)  becomes  small  at  higher  densities,  k, 
becomes  proportional  to  density  with  slope  vcr2.  As  Pp{  1 ) 


L 

FIG.  10.  Pp  (L)  versus  L  for  different  densities. 
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FIG.  11.  Same  as  in  Fig.  2  only  with  fit  parameter  values  determined  from 
their  density-dependent  functions. 

nears  unity  at  low  densities,  k,  becomes  proportional  to  den¬ 
sity  with  slope  vo-\.  This  implies  that  at  higher  pressures,  k, 
does  not  continuously  roll-off  from  an  extrapolation  from  the 
limiting  low-density  results.  There  is  some  evidence  of  this  at 
the  highest  energies  (Fig.  7). 

One  check  of  the  adequacy  of  the  fitted  functions  of  m 
and  k,  as  functions  of  density  is  to  use  the  fitted  values  to 
reproduce  the  decay  curves  in  Fig.  2.  In  Fig.  11,  the  fits  of 
m  and  k,  for  vibrational  relaxation  as  functions  of  density, 
as  shown  in  Figs.  6  and  7,  are  used  to  evaluate  the  LS  form 
of  the  decay  curve  in  comparison  to  the  trajectory -computed 
values.  In  every  case,  except  10  atm,  the  parameter  functions 
of  density  result  in  quite  satisfactory  fits  to  the  trajectory  data. 
For  10  atm,  the  fitted  decay  is  simply  too  slow.  This  is  a 
reflection  of  the  fact  that  in  Fig.  7,  the  largest  relative  fitting 
error  for  kt  occurs  at  the  lowest  density  (corresponding  to  10 
atm).  As  discussed  earlier  (see  Sec.  ID  B),  the  k,  value  in 
Table  I  for  that  density  is  more  representative  of  a  global  fitting 
parameter  than  of  the  initial  decay  rate.  The  actual  initial  rate 
is  noticeably  lower  and  is  represented  by  the  open  circle  in 
Fig.  7.  This  lower  value  of  k,  is  noticeably  closer  to  the  fitted 
value  but  is  still  on  the  high  side.  For  the  10  atm  k,  value 
in  Table  I,  a  satisfactory  relative  error  can  only  be  achieved 
by  using  more  adjustable  parameters  (vcriX  however,  some 
of  the  optimized  values  are  negative  and  therefore  aphysi- 
cal.  The  results  for  k,  at  low  pressures  will  require  further 
studies  and  a  better  understanding  of  the  global  features  of  the 
decay  curves. 

The  optimal  values  of  vcr i  and  vcri  are  0.98  and  0.45, 
respectively.  For  L  >  2,  the  second  number  must  be  multiplied 
by  L  to  get  the  correct  cross  section,  that  is,  vcrL  =  Lmr2.  The 


optimized  values  mean  that  cr2  =  tri  but  crL>2  >  o-\,  one  would 
generally  expect  that  collision  events  involving  multiple  Ar 
atoms  simultaneously  would  more  rapidly  relax  a  molecule 
than  would  a  collision  event  involving  only  one  Ar  atom. 
However,  these  optimized  values  also  mean  that  Lcr\  >  ctl  for 
all  L  >  1.  In  other  words,  the  multi-collider  cross  sections  are 
generally  larger  than  the  single-collider  cross  section  but  not  L 
times  larger.  This  makes  collision  events  involving  more  than 
one  Ar  less  efficient  than  L  separate  collision  events  with  one 
Ar.  This  loss  in  efficiency  makes  k,  undershoot  at  high  densities 
the  extrapolation  from  its  limiting  form  at  small  densities.  This 
undershoot  is  seen  in  Fig.  7  to  occur  at  pressures  starting  at 
~  1 00  atm.  A  plausible  physical  reason  for  this  loss  of  efficiency 
is  that  multiple  Ar  atoms  involved  in  a  collision  expands  the 
number  of  degrees  of  freedom  that  can  absorb  the  collision 
energy,  leading  to  a  “protection”  of  the  target  molecule  by 
“boiling  off”  other  Ar  atoms  involved  in  the  collision.  This 
idea  was  used  by  Schwarzer  et  al. 14  in  an  explanation  of  azu- 
lene/He,Xe  relaxation  measurements. 

The  general  success  of  the  density-dependent  fits  of  m  and 
k\  in  representing  the  computed  vibrational  decay  curves  at  all 
the  higher  densities  motivates  an  extrapolation  of  the  fit  up  to 
solid-state  densities.  The  accuracy  of  this  extrapolation  can  be 
affected  by  the  simple,  statistical  nature  of  the  combinatorial 
model  of  Eq.  (14)  that  renders  all  dynamics  into  two  cross- 
section  parameters  adjusted  to  the  lower  density  results  of  the 
simulations.  Nonetheless,  the  extrapolation  is  instructive  and 
is  shown  in  Fig.  12.  The  solid-state  density  cutoff  depicted 
by  the  dashed  vertical  line  in  the  figure  is  the  packed-hard- 
sphere  limit  discussed  in  Sec.  Ill  B.  The  pressures  needed 
to  condense  Ar  gas  at  300  K  are  very  much  larger  than  the 
pressures  studies  here.  However,  from  the  density  perspective. 


density  (1/A3 ) 


wmmu 

* 

O 

O 


FIG.  12.  The  vibrational  m  and  k,  functions  of  density  in  Figs.  6  and  7 
extended  to  solid-state  densities.  The  dashed  vertical  line  corresponds  to  the 
nominal  solid-state  density.  The  solid  curves  are  the  fits  discussed  in  the  text. 
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the  highest  density  simulated  in  the  present  study  is  already 
25%  of  the  way  to  the  solid-state  cutoff  As  seen  in  Fig.  12, 
both  the  m  and  k,  functions  behave  smoothly  with  density  up 
to  the  solid-state  limit.  The  recovery  of  linear  dependence  of 
the  initial  rate  k,  with  density  at  high  density,  as  discussed  in 
connection  with  Eq.  (13)  and  hinted  at  in  Fig.  7,  is  clearly 
evident  in  the  initial  rate  fit  in  the  approach  to  solid-state 
density.  Because  m  is  substantially  larger  than  unity  over  the 
whole  interval  of  densities  studied,  the  initial  rate  k\  in  Fig.  12  is 
the  maximum  rate  achieved  during  the  decay.  As  more  energy 
drains  out  of  the  molecule,  the  rate  of  subsequent  decay  will 
slow  down. 

For  rotational  relaxation,  as  shown  in  Fig.  8,  the  rate  of 
decay  slows  down  as  rotational  energy  drains  out  of  the  mole¬ 
cule.  However,  unlike  vibrational  decay,  the  k\  for  rotational 
decay  is  essentially  directly  proportional  to  density  for  the 
interval  of  densities  studied.  Rotational  relaxation  occurs  at 
longer  range  than  vibrational  relaxation  because  the  quantum 
of  energy  is  so  much  smaller,  such  that  glancing  collisions  can 
be  effective.  Therefore,  one  would  expect  that  if  simultaneous 
molecular  collisions  with  multiple  Ar  atoms  are  important  for 
vibrational  relaxation,  then  they  are  also  important  for  rota¬ 
tional  relaxation.  In  this  regard,  if  crL  were  the  same  for  all  L, 
then  multi- Ar  collisions  with  nitromethane  could  occur  but  the 
end  result  would  still  be  direct  proportionality  to  density  over 
the  entire  density  interval  examined.  Further  studies  are  needed 
to  clarify  the  differences  between  rotational  and  vibrational 
relaxations. 

The  above  discussion  focuses  on  the  extrapolation  to  the 
high-density  regime  of  the  relaxation  in  the  gas  phase.  Most 
previous  theoretical  studies  of  vibrational  relaxation  are  in  the 
/ow-density  regime  of  isolated  binary  collisions  and  produce 
A E,  defined  earlier  as  the  average  energy  lost  per  collision 
by  the  molecule  with  initial  energy  E.  Unfortunately,  even 
the  lowest  pressure  examined  in  our  study,  10  atm,  is  still  not 
low  enough  to  unambiguously  identify  single  collision  events 
directly  from  the  trajectories.  For  example,  if  one  samples 
the  instantaneous  vibrational  energy  of  nitromethane  every 
tk  -  T+  i  =  50  fs  during  a  typical  10  atm  trajectory  over  an 
extended  period  of  time  (say  100  ps)  and  computes  the  distribu¬ 
tion  of  successive  vibrational  energy  changes,  P(E{t±) 
-E/fk-i)),  this  will  yield  a  distribution  that  has  an  average  very 
near  to  zero,  a  standard  deviation  of  ~0. 1  kcal/mol,  and  a  few 
outliers  that  range  to  as  much  as  ~1  kcal/mol.  A  definition 
of  a  collision  as  one  that  produces  a  change  of  “x”  times 
the  standard  deviation  cannot  be  derived  from  more  rigorous 
theory  and  will  produce  “single”  collision  properties  that  are 
critically  dependent  on  the  value  of  “x.” 

While  the  trajectories  of  our  study  cannot  be  used  to 
unambiguously  identify  single  collision  events,  there  are  ap¬ 
proximate  models  of  Z,  the  number  of  collisions  per  unit 
time,  that  are  based  on  the  cross  sectional  area  of  interaction 
potentials.  Such  models  are  implicitly  derived  in  the  context  of 
isolated  binary  collisions  and  are  directly  proportional  to  the 
bath  gas  density.  Given  a  model  of  Z,  then  in  the  low  pressure 
limit  of  isolated  collisions 


dE  dn  dE 
dt  dt  dn 


(15) 


where  n  is  the  isolated  collision  number  from  the  start  of  the 
trajectories.  From  Eq.  (4)  for  the  LS  functional  form,  it  is  easy 
to  show  that 

dE  t  E(t)\m 

-di-k'm){W))  kiE{0) •  (16) 

Combining  Eq.  (15)  with  the  above  equation  in  the  low- 
pressure  limit  of  isolated  collisions  leads  at  t  —  0  to 

A£  =  -|f(0).  (17) 

Z  is  linearly  dependent  on  pressure  (or  density)  and  in  the  low 
pressure  limit  k\  is  also,  making  their  ratio  and  consequently 
A E  pressure  independent.  The  fit  of  kt  as  a  function  of  density 
follows  from  a  combinatorial  model  that  used  a  hard  sphere 
estimation  of  both  Ar  and  CH3NO2.  For  consistency,  the  well 
known  hard  sphere  model  can  be  used  for  Z 


Z  -  p(RcH3N02  +  Rat) 


2 


8  kT 
n  p 


(18) 


where  u  is  the  reduced  mass  for  Ar  +  CH3NO2.  Using  the  fit  of 
k\  and  values  of  Rch3no2  and  R,\i  derived  previously,  from  Eq. 
(17)  the  low  density  limiting  value  of  AEvn,  is  -0.18  kcal/mol. 

There  are  a  few  single-collision  theoretical  studies  using 
Ar  as  a  bath  gas  where  the  molecule  is  initially  excited  by 
amounts  similar  to  that  studied  here.  The  most  relevant  is 
Lenzer  et  al.4(l  who  studied  rotationally  equilibrated  Q,Hr,  with 
68.6  kcal/mol  initial  vibrational  excitation  in  a  300  K  Ar  bath 
and  found  that  A Evit,=  -0.094  +  0.014  kcal/mol.  The  value 
was  found  to  be  quite  sensitive  to  the  nature  of  the  interaction 
potential,  which  is  not  well  known.  Lendvay  and  Schatz47 
studied  rotationally  equilibrated  SFf,  vibrationally  excited  to 
50  kcal/mol  in  a  300  K  Ar  bath  and  found  that  A Evu,  =  -0.06 
kcal/mol.  Lendvay  et  a l.4''  studied  rotationally  equilibrated 
SO2  vibrationally  excited  to  57  kcal/mol  in  a  300  K  Ar  bath  and 
found  that  A Evu,  =  -0.20  kcal/mol  or  -0.49  kcal/mol  depend¬ 
ing  on  the  nature  of  the  potential.  All  three  of  these  published 
results  involve  a  fixed  initial  vibrational  excitation  and  a  ther¬ 
mal  rotational  distribution,  whereas  the  studies  here  involved  a 
fixed  total  internal  energy  maintained  by  a  correlated  efficient 
microcanonical  sampling  (EMS)  vibrational  and  rotational 
distributions.  Given  the  differences  in  systems  and  initial 
conditions,  our  limiting  single  collision  estimate  of  A F,* 
=  -0.18  kcal/mol  is  not  inconsistent  with  past  studies. 


V.  CONCLUSIONS 

Classical  molecular  dynamics  simulations  were  performed 
to  study  the  relaxation  of  nitromethane  in  an  Ar  bath  as  a  func¬ 
tion  of  pressure  P  =  10, 50, 75, 100, 125, 150, 300,  and  400  atm 
at  300  K.  At  the  highest  pressure,  the  density  of  the  gas  is  ~25% 
of  the  solid-state  density.  The  molecule  was  instantaneously 
excited  by  statistically  distributing  50  kcal/mol  among  the 
internal  degrees  of  freedom.  At  each  pressure,  1000  periodic 
boundary  condition  trajectories  followed  nitromethane  and 
1000  Ar  atoms  explicitly  for  1000  ps,  except  for  10  atm  for 
which  the  simulations  were  extended  to  5000  ps.  The  ensemble- 
averaged  vibrational  and  rotational  energies  were  calculated 
as  functions  of  time.  The  simulations  are  long  enough  for  a 
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substantial  amount  of  the  vibrational  decay,  and  almost  all  of 
the  rotational  decay,  to  occur. 

The  rotational  energy  undergoes  fast  decay  to  ~5%  of  the 
original  energy.  There  follows  a  very  slow  decay,  the  origin 
of  which  is  not  yet  clear  and  for  which  characterization  will 
require  further  study  with  improved  statistics.  The  dominant 
rotational  decay  and  the  vibrational  decay  curves  can  be  satis¬ 
factorily  fit  with  the  Lendvay-Schatz  functional  form,  Eq. 
(4). 19  This  form  has  two  parameters;  one  that  is  sensitive 
to  the  initial  rate  of  decay  and  one  that  is  sensitive  to  the 
curvature  of  the  decay  curve.  All  decay  curves  at  all  pressures 
exhibit  positive  curvature,  which  implies  that  the  rate  of  decay 
decreases  as  energy  drains  out  of  the  molecule  and  is  therefore 
not  describable  by  a  single-exponential  function.  The  initial 
rotational  relaxation  rate  is  nearly  directly  proportional  to 
density  over  the  entire  interval  of  density  studied.  By  contrast, 
the  initial  vibrational  rate  is  not  directly  proportional  to  density 
over  that  same  interval.  Rather,  with  increased  density  the 
initial  rate  falls  away  from  the  extrapolation  of  the  limiting 
low-pressure  proportionality  to  density. 

Fits  to  the  two  Lendvay-Schatz  parameters  were  indepen¬ 
dently  determined  as  functions  of  density  over  the  interval  of 
densities  studied.  For  the  initial  rate,  the  functional  form  arises 
from  a  model  of  the  probability  of  individual  nitromethane 
collision  events  that  involve  multiple  Ar  atoms.  This  model 
is  based  on  the  combinatorial  arithmetic  for  randomly  filling 
space  near  the  nitromethane  molecule  with  Ar  atoms  for 
different  bulk  densities.  The  model  suggests  that  at  the  higher 
densities  studied,  collisions  involving  multiple  Ar  atoms 
become  the  rule  rather  than  the  exception.  Roll-off  of  the  initial 
rate  from  the  low-density  extrapolation  occurs  because  the 
cross  section  for  collision  events  with  L  Ar  atoms  grows  with  L 
more  slowly  than  L  times  the  cross  section  for  collision  events 
with  one  Ar  atom.  The  resulting  density-dependent  function 
of  k,  represents,  reasonably  well,  the  individually  computed 
k,  values,  with  the  poorest  agreement  occurring  at  the  lowest 
density.  Given  the  density-dependent  functions  for  the  two 
parameters  in  the  Fendvay-Schatz  form,  all  the  vibrational 
decay  curves  are  well  represented  except  at  the  lowest  density 
for  which  the  decay  is  somewhat  too  slow.  Extrapolation  of 
the  parameter  functions  gives  a  prediction  of  the  decay  over 
all  densities  corresponding  to  the  gas  phase. 

Further  studies  are  required  to 

•  Characterize  the  residual  rotational  decay; 

•  understand  the  molecular  mechanisms  at  work  for  initial 
vibrational  and  rotational  decays,  in  part,  by  testing  the 
combinatorial  model  for  multiple  collisions;  and 

•  understand  the  molecular  mechanisms  at  work  that  pro¬ 
duce  positive  curvature  for  vibrational  and  rotational 
relaxations. 

The  studies  reported  here  involve  only  one  temperature 
and  only  one  kind  of  initial  molecular  excitation.  A  full  under¬ 
standing  of  the  molecular  mechanisms  involved  will  require 
the  ability  to  predict  and  explain  changes  in  temperature  and 
excitation.  What  the  current  study  has  accomplished  is  the 
development  of  feasible  and  reliable  methods  to  simulate  the 
relaxation  and  a  framework  of  analysis  that  allows  appropriate 
questions  to  be  posed  and  hopefully  answered. 


Additional  studies  of  nitromethane  are  in  the  planning 
stages.  Similar  studies  on  HCb  and  HO  relaxation  in  Ar  are 
in  progress.  Preliminary  results  on  HO2  at  T  -  800  K  indicate 
behavior  qualitatively  similar  to  that  of  nitromethane. 
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